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Analytic  Support  of  the  Phillips  Lab  Whole  Sky  Imager,  1997  -  2001 

Janet  E.  Shields,  Monette  E.  Karr,  Art  R.  Burden, 

Richard  W.  Johnson,  and  Justin  G.  Baker 


1.  Introduction 

The  Atmospheric  Optics  Group  at  the  Marine  Physical  Lab,  Scripps  Institution  of 
Oceanography,  University  of  California  San  Diego,  has  been  developing  Whole  Sky 
Imagers  (WSI)  for  many  years  (Johnson  1989,  1991).  A  Day/Night  WSI  which  had  24- 
hour  a  day  capability  was  initially  developed  in  the  early  1990’s,  and  further  developed 
since  that  time  (Shields  1993,  1998).  One  of  these  instruments,  Day/Night  WSI  Unit  5, 
was  delivered  to  the  Air  Force  Phillips  Lab  in  December  95.  This  work  was  funded 
under  Contract  N00014-93-D-0141  DO  #11  (Shields  1997a). 

A  follow-on  contract,  N00014-96-D-0065-D3  (Shields  1997b)  was  funded  to  support 
additional  documentation  and  to  support  enhancement  of  the  software  provided  with  the 
instrument.  Under  that  contract,  a  new  Operations  Manual  and  Theory  of  Operations 
Manual  were  delivered,  and  the  archival  and  communications  software  was  upgraded  to 
enable  several  new  control  options. 

This  report  documents  the  work  that  was  performed  under  a  second  follow-on  contract, 
Contract  N00014-97-C-0385,  which  was  funded  from  July  97  through  June  01  for  a  total 
of  $  1 03,073.  This  contract  was  intended  to  provide  support  for  additional  upgrades  to  the 
system  software,  algorithms,  and  hardware,  on  an  as-needed  basis.  The  work  was 
supported  by  the  Air  Force  Phillips  Lab,  later  renamed  Air  Force  Research  Lab.  The 
Statement  of  Work  is  given  in  Section  2.  The  initial  work,  to  provide  software  upgrades 
and  hardware  repairs,  is  discussed  in  Section  3.  Much  of  the  funding  increment  was  used 
to  develop  a  first-generation  night  cloud  algorithm,  as  discussed  in  Section  4.  Other 
software  and  hardware  enhancements  performed  during  the  period  when  the  algorithm 
was  under  development  are  discussed  in  Section  5.  The  final  funding  under  this  contract 
was  used  primarily  to  calibrate  the  camera  in  preparation  for  development  of 
transmittance  algorithms.  This  calibration  work  is  discussed  in  Section  6,  which  also 
provides  an  overview  of  ongoing  work  continuing  under  other  more  recent  contracts. 

Although  we  normally  do  not  have  access  to  the  site  and  experiments  where  the 
instrument  is  located,  our  sponsor  described  the  use  of  the  instrument  as  follows:  Air 
Force  Research  Laboratory  employs  the  Whole  Sky  Imager  in  the  role  of  test  support  to 
multiple,  high-value  experiments.  WSI  images  are  displayed  on  web  pages  in  the  test 
control  center  to  increase  situational  awareness.  WSI  data  are  provided  for  both  daytime 
and  nighttime  operations.  During  certain  tests,  the  test  director  will  use  the  Whole  Sky 
Imager  to  either  continue,  or  cancel,  the  test  based  on  the  WSI  cloud  detection  algorithm 
and  resulting  picture.  AFRL  has  supported  continued  algorithm  development  for 
increased  situational  awareness  during  testing  and  collection  of  atmospheric  data  for  post 
processing. 
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2.  Statement  of  Work 


The  Statement  of  Work  from  the  proposal  is  as  follows: 

The  contractor  shall,  unless  otherwise  specified  herein,  supply  the  necessary 
personnel,  facilities,  services,  and  materials  to  accomplish  the  following  tasks  within  a 
two  year  period  following  receipt  of  funding. 

a.  Coordinate  with  the  sponsor  regarding  the  most  appropriate  tasks  and  estimated  costs 
for  development. 

b.  Provide  personnel  trained  in  the  Whole  Sky  Imager  and  its  capabilities  to  address 
these  tasks  (at  MPL)  to  the  limit  of  funding  provided  under  the  contract.  These  tasks  may 
include  analysis,  software  development,  documentation,  minor  hardware  development, 
and  other  tasks  related  to  the  WSI  which  are  mutually  agreed  by  the  sponsor  and  by  MPL 
to  be  appropriate. 

c.  Provide  a  final  written  report  and  interim  verbal  reports  to  the  sponsor  regarding  the 
results  of  the  above  work. 

This  work  was  all  completed  in  a  timely  manner,  i.e.  using  a  schedule  which  best  fit  the 
needs  of  the  sponsor.  The  tasks  were  discussed  in  phone  communications  and  visits 
during  the  period.  The  agreed-upon  tasks  are  discussed  in  the  following  sections. 
Interim  reports  were  provided  in  the  form  of  memos,  phone  discussions,  and  e-mail.  This 
written  report  completes  the  final  requirement.  To  the  best  of  our  knowledge,  the  funded 
work  has  been  completed  to  the  satisfaction  of  the  sponsor. 

3.  Initial  Software  Enhancements  and  Hardware  Repairs,  Jul  97  -  May  98 

The  initial  funding  increment  of  $  10k  was  received  approximately  Jul  97,  and  followed 
by  $20k  in  approximately  Oct  97.  As  agreed  with  the  sponsor,  these  funds  were  used  to 
further  develop  the  software  and  to  repair  the  camera. 

First,  the  program  Viewl6  was  provided  in  order  to  enable  easier  viewing  of  the  WSI 
images.  This  program  also  allowed  the  user  to  save  16  bit  images  as  a  TIFF  image.  As 
documented  in  Technical  Memo  AV98-011t  (Karr  1998a),  an  upgrade  to  the  software 
that  runs  the  WSI,  RunWSI,  was  delivered  in  Feb  98.  This  upgrade  fixed  a  subtle 
problem  between  the  system  clock  and  the  WWV  clock  that  was  occasionally  causing 
problems.  In  addition,  three  “hot  keys”  were  added  for  more  convenient  changes  in  the 
operational  parameters  during  the  program  usage.  Later  upgrades  funded  under  the  later 
funding  increments  are  discussed  in  Section  5. 

Late  in  1997,  the  Photometries  Series  200  16  bit  digital  camera,  which  is  at  the  core  of 
the  Day/Night  WSI,  failed.  We  determined  that  this  was  caused  inside  the  Photometries 
camera  head,  by  corrosion  of  the  metal  surrounding  the  coolant  path.  This  allowed 
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coolant  to  enter  the  chamber  containing  the  camera  CCD  and  the  electronics,  all  of  which 
were  destroyed.  Under  this  contract,  a  new  CCD  chip  was  purchased,  and  CCD  and 
electronics  were  replaced  by  Photometries.  (This  required  separating  the  fiber  optic  taper 
from  the  original  CCD  and  bonding  it  to  the  new  CCD.)  The  full  camera  assembly  (built 
by  MPL  to  hold  the  camera  and  the  optics)  was  then  reassembled,  and  the  camera  was 
recalibrated  and  returned  to  the  field.  This  repair  was  documented  in  Technical  Memo 
AV97-051t  (Shields  1997),  which  also  includes  the  part  numbers  of  the  major 
components. 

The  WSI  camera  housing  was  reinstalled  by  Air  Force  personnel,  however  some 
remaining  problems  resulted  in  a  trip  to  the  site  by  MPL  personnel.  We  were  able  to 
trace  the  problem  to  one  of  the  connectors  installed  by  Photometries  during  the  repair 
(Wertz,  1998). 

4.  Night  Cloud  Algorithm  Development  Work,  May  98  -  Mar  00 

The  next  two  funding  increments,  $20k  received  in  May  98  and  $40k  received  in  May  99, 
were  intended  and  used  primarily  for  the  work  toward  a  night  cloud  algorithm,  with  some 
of  the  funding  used  for  hardware  support  as  needed. 

The  intent  of  the  initial  night  cloud  algorithm  is  to  evaluate  all  regions  of  the  sky  and 
identify  the  presence  or  absence  of  clouds.  (Later  algorithms  which  have  been  developed 
or  are  in  development  address  moonlight  conditions  and  related  issues,  and  should 
provide  optical  depth  or  beam  transmittance.  Concepts  for  the  development  of  a  higher 
spatial  resolution  algorithm  have  also  been  developed.)  This  first  algorithm  works 
essentially  by  evaluating  the  radiometric  signature  of  the  image  in  the  star  field.  We 
initially  wrote  programs  to  analyze  the  character  of  the  image  near  stars,  in  the  presence 
and  absence  of  clouds.  As  reported  informally  in  June  and  July  98  (via  e-mail  and 
phone),  we  evaluated  the  contrast  between  a  star  and  its  background,  and  how  this 
contrast  changed  in  the  presence  of  thin  clouds.  This  initial  study  was  sufficiently 
encouraging  to  continue  the  development. 

4.1  Initial  Star  Mapping 

In  order  to  automate  an  algorithm  based  on  star  detection,  it  is  necessary  to  have  a  means 
of  identifying  the  expected  location  of  stars  in  the  image.  This  requires  the  following: 

a)  a  geometric  calibration  to  convert  from  zenith  and  azimuth  angle  to  x,y  pixel  location. 

b)  an  automated  star  catalog  to  provide  right  ascension  and  declination  of  celestial 
objects 

c)  a  program  to  derive  zenith  and  azimuth  angle  of  these  objects,  given  the  time  and 
location  of  the  WSI,  and  to  convert  them  into  pixel  locations  in  the  image 

Initial  geometric  calibration  equations  were  derived  from  in-house  geometric  calibrations 
as  documented  in  Memo  AV98-041t  (Burden,  1998b).  These  equations  matched  the 
calibration  data  to  within  an  average  error  of  about  .08  degrees  or  about  .2  pixel.  For 
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item  b,  the  US  Naval  Observatory  star  catalog  was  chosen  (Hoffleit  1991).  (This  catalog 
is  also  known  as  the  National  Space  Science  Data  Center  star  catalog.)  This  consists  of 
an  ASCII  file  that  lists  the  right  ascension,  declination,  magnitude,  and  related 
parameters.  A  series  of  programs  were  written  in  IDL1  to  derive  the  positions  of  the  stars 
in  object  space  and  image  space,  as  documented  in  AV98-041t  (Burden,  1998a). 

To  test  the  results,  a  program  was  written  which  extracted  stars  from  the  star  catalog,  and 
plotted  squares  around  the  predicted  location  of  stars  of  magnitude  brighter  than  4.  These 
were  superimposed  on  a  star  image.  We  were  reasonably  pleased  with  the  results, 
however  this  initial  geometric  calibration  was  not  as  accurate  as  we  would  like,  as 
documented  in  Memo  AV98-041t.  The  geometric  calibration  was  refined  as  discussed  in 
the  next  section. 

4.2  Refining  the  Star  Mapping: 

An  image  of  a  star  in  the  WSI  image  plane  is  typically  Gaussian  in  shape,  with  a  width  of 
0.5  pixels.  That  is,  the  point  source  in  object  space  is  convolved  with  a  roughly  Gaussian 
point  spread  function  due  to  the  atmosphere  and  the  optics  in  the  WSI.  Using  this 
information  in  the  image,  one  can  develop  more  precise  techniques  for  determining  the 
star  location.  We  modified  a  technique  developed  by  our  colleague  Dr.  Tim  Tooman 
(verbal  communication),  using  Gaussian  best-fit  routines  to  determine  the  sub-pixel 
location  of  selected  stars.  An  IDL  program  was  written  at  MPL,  GetStarlnfo,  which 
allows  one  to  interactively  select  stars  distributed  over  the  image  and  then  associate  each 
star  with  a  star  in  the  US  Naval  Observatory  catalog.  Using  this  approach,  a  file  is 
created  with  the  stars’  predicted  positions  in  object  space  and  their  measured  positions  in 
image  space. 

The  known  azimuth  and  zenith  angles  of  80  to  100  user-selected  stars,  along  with  the 
decimal  pixel  positions  of  the  stars  in  the  observed  image,  are  then  input  into  a 
Mathematica2  program  Anglefit  written  by  Dr.  Tooman.  This  program  determines  an 
empirical  fit  between  the  object  space  and  image  space  parameters.  Use  of  the  program 
is  somewhat  complex,  and  involves  an  iterative  approach,  including  evaluating  the 
residuals.  It  is  slightly  subjective,  in  that  it  allows  the  user  to  estimate  the  center  position 
(actually  the  position  on  the  image  corresponding  to  the  zenith),  and  reruns  are  made 
until  the  user  is  satisfied  with  the  residuals.  Results  are  typically  to  within  an  average 
error  of  0.5  pixels  or  better.  The  resulting  geometric  equations  are  documented  in 
Technical  Memo  AV99-043t  (Burden  1999a). 


1  IDL,  which  stands  for  "Interactive  Data  Visualization",  is 
produced  by  Research  Systems,  Inc.,  which  was  recently 
purchased  by  Kodak.  Their  address  is  4990  Pearl  East 
Circle,  Boulder,  CO  80301. 

2Mathematica  is  produced  by  Wolfram  Research,  Inc.,  whose 
address  is  100  Trade  Center  Drive,  Champaign,  IL  61820. 
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4.3  Development  of  the  Algorithm 


The  next  step  in  the  algorithm  development  was  to  develop  a  program  to  detect  stars  in  a 
test  image.  The  program  uses  the  US  Naval  Observatory  star  catalog  to  determine  the 
stars  within  a  selected  magnitude  range,  and  predicts  their  anticipated  location  in  a  test 
image  (in  image  space).  The  program  then  finds  the  largest  pixel  value  within  a  two- 
pixel  radius  of  the  predicted  location.  Using  a  Levenberg-Marquardt  best-fit  routine3,  the 
program  attempts  to  fit  a  two-dimensional  Gaussian  curve  to  the  data  (and  in  the  process 
determines  the  true  pixel-center  of  the  peak  in  image  space).  If  the  fitting  routine  gives  a 
reasonable  result,  the  program  counts  this  as  a  successful  detection  of  the  star  in  the 
image.  Using  this  technique,  we  were  able  to  detect  most  of  the  stars  from  magnitude  1 
to  magnitude  7  in  a  sample  clear  image.  For  example,  for  zenith  angles  30  degrees  or 
less,  the  program  detected  100%  of  the  33  stars  of  magnitude  0  to  4,  and  88%  of  the  99 
stars  of  magnitude  0  to  7.  For  zenith  angles  of  70  degrees  or  less,  the  program  detected 
78%  of  the  stars  of  magnitude  0  to  7.  Many  of  the  undetected  stars  were  automatically 
rejected  by  the  algorithm,  because  it  rejects  any  stars  which  are  near  other  stars.  A  test  of 
random  locations  in  an  overcast  image  resulted  in  a  5%  false  detection  rate. 

On  the  basis  of  the  above  results,  a  first  generation  cloud  algorithm  was  generated  as 
documented  in  Technical  Memo  AV99-044t  (Burden,  1999b),  which  is  included  in 
Attachment  1.  This  program  divides  the  sky  out  to  a  70  degree  zenith  angle  into 
segments  which  are  approximately  5  degrees  in  zenith  angle  by  15  degrees  azimuth 
angle.  (Near  the  center  fewer  azimuthal  segments  are  made.)  This  results  in  237  regions. 
Then  within  each  region  the  program  determines  whether  there  are  any  stars  from  the  star 
catalog,  and  if  so,  determines  whether  they  can  be  detected  simply  by  the  presence  of 
their  PSF.  The  program  marks  the  region  grey  if  there  are  no  candidate  stars  in  the  star 
catalog,  blue  more  than  50%  of  the  available  stars  in  the  region  are  detected,  and  red  if 
fewer  than  50%  of  the  available  stars  are  detected. 

This  approach  worked  reasonably  well,  as  shown  in  the  memo.  On  the  clear  image,  90% 
of  the  regions  had  candidate  stars,  and  83%  of  the  regions  with  candidate  stars  were 
identified  as  clear.  The  overcast  image  had  a  thin  overcast;  Orion  could  be  seen  through 
the  clouds,  but  most  of  the  dimmer  stars  could  not  be  seen  in  the  image.  In  this  image, 
88%  of  the  regions  had  candidate  stars,  and  3%  of  the  regions  with  candidate  stars  were 
identified  as  clear.  The  scattered  cloud  case  was  in  between,  with  66%  identified  as 
clear. 
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For  the  LM  best  fit  technique,  we  used  a  routine,  MPFIT  IDL 
software,  written  by  Craig  Markwardt  at  NASA.  It  was 
obtained  from  Craig  B.  Markwardt,  NASA/ GSFC  Code  662, 
Greenbelt,  MD  20770,  craigm@ lheamail . gs f c . nasa . gov, 
http: / / cow.physics.wisc.edu/-craigm/idl/idl.html .  For  a 
further  reference,  consider 

Press,  W.H.,  Et .  A1 . ,  1992,  Numerical  Recipes  In  Fortran, 
Cambridge  University  Press,  p.  678-682. 


5 


At  this  point,  there  were  several  productive  directions  we  might  go  in.  The  decision  was 
made  to  convert  the  program  in  its  current  (at  that  time)  state  to  C  code,  in  the  hopes  of 
fielding  it  with  the  instrument  at  a  later  time.  Other  potential  development  are  discussed 
in  Section  4.5. 

4.4.  Conversion  of  the  Starlight  Algorithm  to  C 

The  last  step  of  the  algorithm  development  funded  under  this  contract  was  the  conversion 
of  the  code  from  the  interactive  IDL  code  developed  as  described  in  Section  4.3  into  a 
more  automated  C  code.  This  version  made  more  use  of  input  files,  as  opposed  to  user 
input,  and  was  considerably  faster  (on  the  order  of  seconds  as  opposed  to  minutes  to 
process  an  image).  This  version  of  the  code  is  documented  in  Memo  AV99-054t  (Burden 
1999e). 

4.5  Further  Development  of  the  Starlight  Algorithm 

Although  the  sponsor  did  not  have  sufficient  funds  to  complete  the  algorithm  in  a  version 
which  could  be  fielded  with  the  instrument,  this  work  has  been  completed  under  funding 
from  another  sponsor.  During  this  later  work,  the  algorithm  was  also  improved 
considerably  in  sophistication,  to  enable  it  to  handle  moonlight  and,  to  some  extent, 
nearby  city  lights. 

At  present  the  algorithm  is  based  on  uncalibrated  raw  data.  There  are  several 
enhancements  which  should  be  feasible  if  using  calibrated  data,  and  are  invarious  stages 
of  development.  First,  it  should  be  quite  feasible  to  provide  estimates  of  earth-to-space 
beam  transmittance  or  optical  depth.  By  evaluating  the  calibrated  signal  from  a  large 
number  of  stars,  and  comparing  it  with  the  return  expected  from  the  stars  in  the  absence 
of  clouds,  we  feel  that  we  can  develop  a  means  for  providing  a  map  of  the  variance  in 
absolute  earth-to-space  beam  transmittance.  This  work  has  been  funded  under  a  separate 
contract. 

In  addition,  better  handling  of  moonlight  and  of  city  lights  should  be  possible.  Our 
understanding  from  our  sponsors  who  are  using  the  automated  algorithm  is  that  they  are 
quite  pleased  with  the  results,  but  we  have  not  yet  done  extensive  analysis  of  the  results 
and  the  potential  for  improvement. 

Also,  we  have  developed  approaches  for  improving  the  resolution  of  the  night  algorithm 
to  potentially  provide  pixel-by-pixel  resolution.  While  we  have  not  been  funded  to 
address  this  work,  we  feel  that  it  is  a  definite  possibile  future  application  for  the  WSI 
night-time  data.  An  additional  useful  development  would  be  the  application  to  sunrise 
and  sunset  time  regimes. 

5.  Additional  Work  Performed  during  the  May  98  -  Mar  00  Period 

Although  the  funding  during  this  period  was  primarily  directed  toward  night  algorithm 
development,  as  discussed  in  Section  4,  additional  work  to  support  the  ongoing 
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operational  use  of  the  WSI  was  completed  as  needed.  Another  software  upgrade  was 
delivered  in  Oct  98,  as  documented  in  Technical  Memo  AV98-053t  and  55t  (Karr  1998b 
and  c).  This  upgrade  included  addition  of  a  time  sync  operation  that  allows  the  display 
computer  to  synchronize  with  the  collection  computer.  Also,  a  SCRAMnet  life  check 
feature  was  added  to  provide  a  way  to  check  if  the  SCRAMnet  collection  is  responding. 
Additional  features  were  added  as  documented  in  the  tech  memo.  We  also  evaluated 
changes  required  for  Y2K  handling  and  made  the  necessary  changes.  An  additional 
minor  upgrade  is  documented  in  Memo  AV99-05 1 1  (Karr,  1 999). 

The  sponsor  decided  that  they  would  like  to  use  the  new  geometric  calibration  results  to 
improve  the  accuracy  of  the  angle  assignment  to  the  pixels  in  the  field  images.  Since 
Unit  5  is  at  a  classified  site,  we  are  restricted  from  viewing  and  processing  images.  For 
this  reason,  we  developed  the  geometric  calibration  procedures  into  a  series  of  steps  that 
could  be  performed  by  the  sponsor,  and  wrote  and  delivered  two  tutorial  memos 
documenting  these  processes.  These  memos  are  Technical  Memoranda  AV99-049t  and 
050t  (Burden,  1999c  and  d). 

A  new  occultor  trolley  motor  was  purchased  during  this  interval,  and  installed  at  the  site 
in  August  99,  as  documented  in  Technical  Memo  AV99-046t  (Wertz  and  Shields  1999). 

6.  Preparation  for  Beam  Transmittance  Development  Work,  Mar  00  -  Jun  01 

The  final  funds  transfer  of  $13,073  was  received  approximately  Mar  00.  At  that  time,  the 
request  was  that  these  funds  be  saved  for  camera  repair,  in  case  any  repair  was  required. 
No  repairs  were  required  initially.  In  the  meantime,  considerable  development  of  the 
night  algorithm  occurred  under  funding  from  other  sponsors,  with  the  result  that  the  night 
algorithm  was  fielded  on  other  Day/Night  WSI  units.  Following  meetings  with  the 
sponsor  in  October  00,  the  decision  was  made  to  fund  the  enhancement  of  the  night 
algorithm  to  determine  earth-to-space  beam  transmittance.  This  task  was  funded  under  a 
follow-on  contract,  N-00014-01-D-0043  DO  #5.  Development  of  the  beam  transmittance 
algorithm  would  require  calibration  of  the  camera,  so  the  remaining  funds  under  the 
existing  contract  were  used  to  acquire  the  radiometric  calibration  data  for  the  camera. 
The  camera  was  also  repaired  at  this  time,  as  documented  in  Memo  AV01-040t  (Baker, 
2000).  Initial  work  on  the  calibration  processing  was  also  begun,  and  although  the 
calibration  processing  was  completed  under  the  follow-on  contract,  the  calibration 
memos  AV01-071  -  079  are  referenced  here.  Work  on  the  beam  transmittance 
development  is  in  progress  at  this  time,  under  the  follow-on  contract. 

7.  Summary 

The  raw  WSI  night  imagery  provides  extremely  useful  images,  which  are  easily 
interpreted  by  human  users.  As  may  be  seen  from  the  sample  images  in  the  attachment, 
both  clouds  and  stars  can  be  readily  seen,  and  the  general  character  of  the  clouds  is 
readily  observed.  To  the  best  of  our  knowledge,  there  is  no  other  instrument  in  existence, 
which  has  this  combination  of  features,  allowing  human  interpretation  of  the  cloud  scene 
during  night-time,  daylight,  and  in  sunrise  and  sunset  time  regimes.  In  addition,  because 
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the  data  have  excellent  dynamic  range  and  low  noise,  and  can  be  calibrated  for  absolute 
radiance,  they  are  very  useful  for  further  mathematical  evaluation.  The  starlight 
algorithm  is  an  example  of  a  first-generation  algorithm  which  starts  to  utilize  this 
outstanding  data. 

During  the  period  from  July  97  through  June  01,  Contract  N00014-97-C-0385  provided 
funding  to  support  the  use  of  the  Day/Night  Whole  Sky  Imager  by  Air  Force  Phillips  Lab 
(later  renamed  Air  Force  Research  Lab).  The  sponsors  were  very  happy  with  the 
instrument,  reporting  to  us  that  their  program  “has  been  very  happy  with  the  WSI  system 
we  purchased  from  MPL  a  few  years  back.  It  is  an  excellent  source  of  data  and  is  a  very 
reliable  system.”  (personal  e-mail).  During  this  time,  repairs  were  performed  as  needed, 
and  software  upgrades  were  delivered.  A  means  of  providing  a  very  accurate  assignment 
of  the  angular  location  (in  object  space)  of  each  pixel  (in  image  space)  was  developed 
and  delivered.  A  first  generation  night  algorithm,  for  identifying  the  presence  of  clouds 
in  a  starlight  image,  was  developed.  Preliminary  work  toward  extending  these  techniques 
toward  determination  of  earth-to-space  beam  transmittance  was  begun,  and  is  continuing 
under  follow-on  funding.  This  work  has  enabled  the  Air  Force  to  continue  many 
productive  years  of  using  the  WSI  in  a  variety  of  test-support  scenarios. 
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MARINE  PHYSICAL  LABORATORY,  0701 
of  the  Scripps  Institution  of  Oceanography 
San  Diego,  California  92152-6400 


AV99-044t 

17  August  1999 


This  memo  describes  preliminary  work  toward  developing  a  nighttime,  starlight-based 
cloud  algorithm  to  be  included  in  the  day-night  WSI  system.  Related  information  regarding 
automatic  detection  of  stars  in  nighttime  imagery  can  be  found  in  Tech.  Memos  AV98-040t, 
AV98-042t,  and  AV99-043t. 

Background 

The  day/night  WSI  has  test  software  capable  of  detecting  clouds  at  night  under  full-moon 
conditions,  but  no  software  is  available  for  detecting  clouds  at  night  under  moonless  conditions. 
We  are  currently  working  on  development  of  a  moonless  night  cloud-detection  algorithm  for  use 
with  the  WSI.  Here,  we  discuss  some  initial  findings  and  present  a  preliminary  cloud-detection 
algorithm  written  in  IDL.  The  current  capabilities  of  the  algorithm  are  discussed,  as  well  as 
issues  that  we  feel  need  further  development. 

The  approach  of  the  initial  study  has  been  to  rely  on  the  difference  in  the  optical 
characteristics  of  a  star  under  cloud-free  conditions  compared  to  cloudy  conditions.  When  an 
image  is  formed  by  an  optical  system,  the  image  is  normally  slightly  defocused.  Mathematically, 
the  image  is  convolved  with  the  point-spread  function  (PSF)  of  the  optics.  For  example,  the 
image  of  a  point  source,  such  as  a  star,  may  be  spread  out  over  several  pixels.  The  PSF  for  the 
WSI  optical  system  can  be  modeled  by  a  Gaussian  distribution.  Using  images  of  stars,  we  have 
found  that  the  width  of  the  PSF,  which  is  equal  to  the  standard  deviation  of  the  Gaussian 
distribution,  is  typically  0.3-0.8  pixels.  It  does  not  appear  to  vary  systematically  over  the  image. 

An  image  of  a  star  centered  on  a  pixel  is  characterized  by  a  peak  signal  at  the  center  pixel 
surrounded  by  pixels  exhibiting  signals  that  decrease  to  the  background  value  further  from  the 
peak.  A  star  with  a  point-spread  function  of  0.5  is  characterized  by  a  signal  that  is  reduced  to 
background  within  three  pixels  of  the  center.  If  the  same  star  is  obscured  by  opaque  cloud,  most 
of  the  signal  from  the  star  will  not  be  transmitted  through  the  cloud  to  the  camera.  The  signal 
from  the  cloud  will  not  exhibit  Gaussian  behavior,  and  will  most  likely  exhibit  only  small 
variations  in  signal.  If  the  cloud  covering  the  star  is  semi-transparent,  a  reduced  signal  from  the 
star  may  still  show  through.  The  background  signal  may  be  increased  by  the  presence  of  the  thin 
cloud. 


Technical  Memorandum 

To:  Atmospheric  Optics  Group 
From:  A.R.  Burden  and  J.E.  Shields 

Subject:  Preliminary  Results  of  Starlight-Based  WSI  Cloud  Algorithm 


Using  a  star  catalog  to  determine  the  locations  of  bright  stars  within  a  WSI  image,  we  can 
search  for  the  stars  using  tests  based  on  the  differences  in  optical  characteristics  just  described. 
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The  presence  of  a  Gaussian-like  peak  at  the  expected  star  location  can  be  used  as  an  indicator  of 
cloud-free  conditions,  while  the  absence  of  the  peak  can  be  used  as  an  indicator  of  cloud 
obscuration.  This  basic  feature,  along  with  other  tests  based  on  differences  in  signal  for  stars  and 
clouds,  is  incorporated  into  the  first-generation  cloud-detection  algorithm  described  below.  This 
scheme  represents  one  possible  method  for  identifying  clouds  in  WSI  imagery.  Other  possible 
schemes  may  be  considered  as  funding  permits. 

The  algorithm  described  here  has  not  been  written  in  a  form  that  is  ideal  for  real-time  use 
with  the  WSI.  Rather,  it  was  written  in  IDL  to  facilitate  the  development  of  such  an  algorithm 
that  could  later  be  written  in  C  and  incorporated  into  the  field  software.  The  basic  structure  of 
the  current  algorithm  can  be  divided  into  two  parts:  subroutines  that  obtain  information  about  the 
catalog  stars  and  the  image  pixels  that  correspond  with  these  stars,  and  subroutines  that  use  this 
information  in  tests  that  differentiate  cloud-free  and  cloudy  regions. 

The  first  group  of  subroutines  is  called  by  the  main  routine  GetStarData.  The  IDL 
version  of  GetStarData  can  be  very  time  consuming  to  run,  as  it  stores  data  for  thousands  of  stars 
in  an  image  with  magnitudes  in  the  range  0-7.  The  information  from  GetStarData  is  saved  in  a 
file  for  each  image,  and  available  for  future  analysis  using  the  second  part  of  the  algorithm.  At 
this  time,  the  cloud  decision  part  of  the  algorithm  is  found  in  a  routine  called  StarAlgorithml. 
As  new  methods  of  cloud  detection  are  developed,  they  can  be  incorporated  into  the  algorithm 
and  applied  to  GetStarData  output  files  without  having  to  rerun  the  entire  set  of  procedures. 

Data 


In  order  to  develop  and  test  the  new  nighttime  algorithm,  images  have  been  obtained 
from  an  archive  of  data  collected  with  the  WSI  located  at  the  SGP  site  operated  by  ARM.  These 
images  were  obtained  by  contacting  the  ARM  archive  and  requesting  data  from  several  nights 
without  moonlight  during  the  years  1997-1999.  The  data  were  downloaded  from  the  ARM 
archive  web  site  and  extracted  from  compressed  archive  files.  Each  image  file  was  then 
calibrated  using  the  IDL  software  Wsi_SGP_Data_Ingest,  written  by  Tim  Tooman  of  Sandia 
National  Laboratory,  using  the  radiometric  calibration  equations  developed  at  MPL.  The 
resulting  files  are  in  the  NetCDF  data  format.  Files  written  in  the  NetCDF  format  can  be  read 
using  standard  IDL  routines.  The  MPL-written  routines  GetCDFData  and  GetCDFDate  are  used 
to  extract  the  image  and  header  information  for  use  in  the  current  study. 


GetStarData 

Locating  Star  Positions 

In  order  to  identify  a  star  in  WSI  imagery,  the  theoretical  location  of  the  star  and  a 
transformation  from  the  theoretical  coordinates  to  image  coordinates  are  required.  The  first  item 
is  easily  obtained  from  one  of  many  bright  star  catalogs  available.  Star  catalogs,  such  as  the 
National  Space  Science  Data  Center  (NSSDC)  bright  star  catalog  used  here,  contain  vital 
information  about  known  stars,  including  right  ascension,  declination,  magnitude,  and  reference 
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numbers.  The  NSSDC  catalog  used  here  is  a  formatted  ASCII  file  containing  over  9100  stars  of 
up  to  magnitude  8  (magnitude  decreases  with  brightness). 

The  geometric  calibration  routine,  Rad2Az,  described  in  AV98-040t  is  used  to  convert 
celestial  coordinates  (right  ascension,  declination)  for  a  star  to  spherical  coordinates  (azimuth, 
zenith)  when  the  time  of  observation  is  known.  Although  coarse  equations  to  derive  pixel 
locations  from  spherical  coordinates  are  given  in  AV98-040t,  we  now  use  a  more  precise 
geometric  calibration  based  on  an  empirical  best  fit  that  is  a  function  of  azimuth  and  zenith 
angles  to  obtain  image  pixel  coordinates.  This  procedure,  described  in  AV99-043t,  is  capable  of 
locating  stars  in  the  catalog  to  within  half  a  pixel.  The  conversion  routine  that  uses  the  empirical 
fit  equations  is  called  GetStarXY. 

The  conversion  equations  obtained  with  the  fitting  procedure  must  be  reworked  for 
different  WSI  cameras  and  each  time  a  WSI  camera  is  moved.  In  the  case  of  SGP  data,  the 
camera  was  moved  on  a  couple  of  occasions  during  the  period  of  study,  December  18,  1997  to 
January  17,  1999.  Therefore,  three  sets  of  calibration  equations  were  developed,  and  the 
appropriate  set  is  used  according  the  date  given  in  the  filename  of  each  image. 

Identifying  Stars  in  Cloud-Free  Images 

Once  the  position  of  a  star  in  image  pixel  coordinates  has  been  determined,  the  next  step 
is  to  determine  whether  a  star-like  peak  exists  in  the  image  at  that  point.  For  this  study,  we  use 
an  iterative  method  that  searches  the  region  surrounding  the  calculated  star  center  for  a  peak. 
First,  the  maximum  signal  for  the  pixels  directly  next  to  and  including  the  calculated  star  center 
is  determined.  If  the  pixel  with  maximum  signal  is  the  calculated  star  center,  the  search  is 
stopped.  Otherwise,  the  pixel  with  maximum  signal  is  labeled  as  the  new  center  pixel  and  the 
process  is  repeated.  If  no  star  center  can  be  found  within  a  certain  pixel  distance  of  the 
calculated  star  center,  the  search  is  considered  a  failure.  For  this  study,  a  maximum  distance  of 
three  pixels  was  used.  Because  crowded  star  fields  may  contain  multiple  star  peaks  in  the  same 
region,  a  second  attempt,  more  sophisticated  attempt,  is  made  to  find  a  different  peak  near  the 
calculated  star  center  when  the  first  attempt  fails.  The  IDL  routine  that  searches  for  potential  star 
peaks  is  called  FindCenterPixel. 

Fitting  Potential  Stars  to  2-Dimensional  Gaussian  Distribution 

If  a  potential  star  peak  is  found  near  the  calculated  star  position,  a  best  fit  calculation  is 
made  using  the  peak  pixel  and  the  surrounding  pixels.  The  best  fit  subroutine  uses  the 
Levenberg-Marquardt  method  (LMM),  which  is  a  robust  technique  commonly  used  for  solving 
non-linear  least-squares  problems  (see  Numerical  Recipes,  Press,  et  al,  for  instance).  Additional 
information  on  the  method  will  be  documented  in  a  future  memo  describing  software  for 
focusing  the  WSI  camera.  The  array  of  pixel  values  that  represents  the  star  are  fitted  with  a  two- 
dimensional  Gaussian  function  with  six  parameters,  as  shown  in  Equation  1.  An  additional 
parameter,  the  tilt,  which  can  be  used  to  determine  the  orientation  of  an  asymmetrical  Gaussian 
about  the  axes,  has  been  left  out  for  now.  While  the  asymmetric  nature  of  a  star  distribution  is 
more  realistic,  model  results  to  date  have  not  been  reasonable. 
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with  the  following  parameters: 

A0  =  Background  value 
A,  =  Scaling  factor 

A2  =  Width  of  Gaussian  in  the  x-direction  (point  spread) 
A3  =  Width  of  Gaussian  in  the  y-direction  (point  spread) 
A4  =  Center  in  x-direction 
A5  =  Center  in  y-direction 


LMM  finds  the  set  of  parameters  that  best  fits  the  array  of  pixel  values.  The  fit  is  "best" 
in  the  least-squares  sense;  that  is,  the  sum  of  the  weighted,  squared  differences  between  the 
model  and  data  is  minimized.  While  calibrated  data  were  used  in  this  case,  the  routines  work 
similarly  for  uncalibrated  WSI  images.  The  IDL  subroutines  that  perform  the  best  fit  calculation 
are  called  MPFit  and  MPFitFun.  These  subroutines  were  converted  from  Fortran  routines 
included  in  the  MinPack  statistical  package  by  Craig  B.  Markwardt  of  NASA.  Additional 
information  can  be  obtained  at  his  website:  http://cow.physics.wisc.edu/~craigm/idl. 

As  described  in  the  Validation  section,  the  LMM  has  several  optional  inputs  that  may 
affect  the  results:  an  array  may  be  passed  to  MPFitFun  that  contains  weights  for  each  pixel  value; 
the  six  parameters  can  be  given  constraints  on  their  values;  and  a  ceiling  can  be  put  on  the 
number  of  iterations  performed.  Additionally,  multiple  best  fit  calculations  may  be  required 
when  the  original,  5-  by  5-pixel  array  does  not  provide  a  good  fit.  The  1 -sigma  errors  produced 
by  the  covariance  matrix  of  the  best  fit  are  used  as  an  indicator  of  good  fit,  and  may  be  high 
when  the  array  size  is  not  large  enough.  When  a  parameter  error  is  higher  than  a  certain 
threshold,  the  fit  is  considered  suspect.  In  this  case,  the  array  width  is  increased  by  2  pixels,  and 
the  fit  is  performed  again.  A  maximum  array  size  of  1 1  by  1 1  pixels  is  used.  While  somewhat 
arbitrary,  this  method  has  shown  to  give  better  results  than  using  either  small  arrays  or  large 
arrays  alone.  Some  stars  require  additional  surrounding  pixels  to  be  distinguished  from  the 
background,  while  others  lie  in  crowded  star  fields  where  too  large  an  array  prohibits  a  good  fit 
with  the  Gaussian  model. 

For  each  best  fit  calculation,  important  information  is  stored  regarding  the  results.  As 
stated  in  the  previous  paragraph,  when  a  best  fit  calculation  results  in  large  error  values,  a  new 
calculation  is  made  with  a  larger  input  array.  Other  factors  that  merit  a  new  calculation  include 
iterations  that  exceed  a  specified  maximum  and  star  point-spread  function  widths  that  exceed 
specified  limits.  For  this  study,  calculations  exceeding  300  iterations  or  point-spread  function 
widths  outside  the  range  0-2  were  performed  again  with  larger  arrays.  Most  best  fit  calculations 
converge  within  a  few  iterations,  especially  for  brighter  stars.  If  no  array  width  between  5  and 
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1 1  pixels  results  in  an  acceptable  outcome,  the  star  in  question  is  flagged  with  an  appropriate 
code  (Code  1  if  maximum  iterations  or  point-spread  limits  exceeded;  Code  4  if  large  errors) 
identifying  the  cause  of  the  failure.  Additional  flags  are  set  if  the  region  near  the  calculated  star 
center  is  potentially  crowded  with  stars  (Code  2:  multiple  peaks  nearby),  or  if  the  star  peak 
coincides  with  a  previously  identified  star  (Code  3:  stars  with  identical  center).  All  other 
outcomes  are  classified  as  a  success  (Code  0).  The  list  of  MPL-written  subroutines  used  in  the 
best-fit  calculation  includes  CheckForPeak,  StarPeak,  RunMPFit,  and  Gauss2D_Mod. 

Statistics  and  Data  Storage 

The  next  step  after  the  best  fit  has  been  determined  is  to  calculate  certain  statistics  for  the 
star,  or  lack  thereof.  The  contrast  is  calculated  by  finding  the  integral  of  the  two-dimensional 
Gaussian  model  of  the  star  and  dividing  the  result  by  the  background  value  for  the  star.  The 
background  is  just  the  first  parameter  resulting  from  the  best  fit.  The  peak  is  equal  to  the  sum  of 
the  background  and  scaling  factor  (i.e.,  Aq+A^.  The  full-width-half-maximum  (FWHM)  is 
defined  for  a  Gaussian  distribution  to  be  •v/8Ln(2)a  ,  or  2.36af.  Since  in  this  case  the  distribution 
is  two-dimensional,  we  take  the  mean  of  the  x  and  y  FWHM  values. 

For  cases  where  no  satisfactory  fit  was  found,  the  statistics  are  calculated  differently. 
The  background  is  the  median  of  the  9-  by  9-pixel  array  surrounding  the  calculated  star  center. 
The  peak  is  the  mean  of  the  3-  by  3-pixel  array  surrounding  the  calculated  star  center.  The 
contrast  is  the  peak  minus  the  background  divided  by  the  background.  The  FWHM  is  set  equal 
to  0.  The  statistical  calculation  routines  are  StoreStats  and  Gauss2Dintegral. 

Currently,  for  each  image  analyzed,  information  is  recorded  for  all  stars  in  the  NSSDC 
catalog  with  magnitudes  in  the  range  0-7  and  zenith  angles  less  than  70°.  The  following 
information  is  stored  for  each  star:  HR  and  SAO  catalog  numbers,  right  ascension,  declination, 
magnitude,  multiple  star  designation,  zenith,  azimuth,  expected  x-pixel,  expected  y-pixel,  image 
x-pixel,  image  y-pixel,  quality  flag,  background,  peak,  contrast,  the  x  and  y  point-spread 
functions,  and  a  image  cell  number  described  in  the  next  section. 


StarAlgorithml 

Below  is  a  description  of  the  preliminary  version  of  the  nighttime,  starlight-based  cloud 
algorithm  based  on  initial  study  results.  While  simple  in  its  application,  it  demonstrates  that  the 
basic  concepts  described  in  the  previous  sections  can  yield  reasonable  results.  The  approach 
used  here  is  to  define  a  number  of  regions,  or  ‘image  cells’,  that  are  approximately  the  same  size 
and  oriented  around  the  image  center,  as  shown  in  Figure  1.  The  grid  is  made  up  of  a  total  of 
237  image  cells  covering  zenith  angles  up  to  70°.  Each  image  cell  is  examined  to  determine  the 
number  of  potential  stars  that  exist  there.  If  at  least  half  the  potential  stars  are  identified  as 
exhibiting  star-like  qualities,  the  entire  image  cell  is  classified  as  cloud-free. 


*  The  Full  Width  at  Half  Maximum  is  the  width  of  the  Gaussian  distribution  at  the  half-power  point,  which  is  equal 
to  0.5  for  a  normal  distribution.  The  value  shown  here  is  obtained  by  solving  the  one-dimensional  Gaussian  for  x 
and  multiplying  the  result  by  2. 
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Figure  1 .  WSI  Image  Overlayed  with  Cloud  Decision  Grid 


Parameters  that  are  considered  important  for  identifying  a  star  include:  the  contrast  of  the 
star  peak  with  the  sky  background;  the  width  of  the  point  spread  functions  associated  with  the 
Gaussian  model  of  the  star;  and  the  value  of  the  quality  flag  based  on  the  best  fit  determination 
for  the  star.  Additionally,  the  contrast  of  a  star  as  well  as  the  ability  of  the  best  fit  algorithm  to 
successfully  converge  is  diminished  for  higher  magnitude  stars.  Therefore,  the  quality  flag, 
contrast,  and  point-spread  functions  of  low  magnitude  stars  within  an  image  cell  are  examined 
first.  If  none  exist,  then  all  other  stars  are  examined.  The  magnitude  threshold  used  here  is  5. 

Before  the  clear  and  cloudy  regions  in  the  image  are  identified,  all  stars  that  were  flagged 
as  being  stacked  on  other,  lower  magnitude  stars  (Code  3),  are  removed  from  the  dataset.  Each 
image  cell  is  then  categorized  as  cloud-free  or  cloudy.  The  first  step  of  the  decision  process  is  to 
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identify  any  star  in  the  image  cell  with  Code  1,  defined  here  as  a  best  fit  failure,  as  cloudy.  Next, 
all  stars  in  the  image  cell  with  a  contrast  >0.18  and  x-  and  y-point  spread  function  widths  in  the 
range  0.3-0.8  are  identified  as  cloud-free.  All  other  stars  in  the  image  cell  are  classified  as 
cloudy.  If  the  number  of  cloud-free  stars  is  at  least  half  the  total  number  of  low  magnitude  stars 
in  the  cell,  the  entire  cell  is  identified  as  cloud-free.  If  no  low  magnitude  stars  are  found  in  the 
image  cell,  the  same  process  is  done  for  the  high  magnitude  stars.  At  this  time,  the  additional 
information  from  Codes  2  and  4  is  not  used. 

The  thresholds  used  in  the  algorithm  were  obtained  by  examining  the  point-spread 
function  widths  of  stars  in  completely  cloud-free  images  and  by  comparing  the  contrast  of  stars 
in  completely  cloud-free  images  and  completely  overcast  images.  Figure  2  shows  the 
distribution  of  point-spread  functions  for  stars  identified  in  several  cloud-free  images.  The  data 
were  separated  into  low  and  high  magnitude  stars,  and  stars  flagged  as  having  poor  best  fit 
convergence  were  removed  from  the  data.  Each  distribution  was  fitted  with  a  Gaussian 
distribution  to  estimate  the  mean  and  standard  deviation.  In  both  cases,  the  point-spread  function 
widths  were  centered  near  0.48,  with  most  cases  being  in  the  range  0.3-0.8.  The  high  magnitude 
distribution  exhibited  a  slightly  lower  mean,  but  skewed  toward  higher  point-spread  function 
widths.  This  demonstrates  the  difficulty  in  performing  best  fit  calculations  on  high  magnitude 
stars. 


Figure  3  shows  contrast  distributions  for  stars  under  cloud-free  conditions  and  stars  under 
cloudy  conditions.  Again,  stars  were  divided  into  low  and  high  magnitude  datasets.  In  order  to 
highlight  the  differences  in  the  distributions,  only  contrasts  up  to  a  value  of  5  were  plotted.  In 
the  cloud-free,  low  magnitude  case,  most  of  the  values  were  in  the  range  0.25-2.5.  A  few  stars 
had  contrasts  as  high  as  50  or  more.  The  contrasts  for  the  cloud-free,  high  magnitude  case  were 
much  lower,  with  most  of  the  values  lying  in  the  range  0.05-1.0.  Under  completely  overcast 
conditions,  the  contrasts  rarely  increase  above  0.2,  and  these  cases  may  be  due  to  optical 
aberrations,  bright  stars  being  viewed  through  clouds  with  lower  optical  depth,  or  natural 
variations  in  cloud  radiance.  Since  some  crossover  exists  between  the  cloud-free  and  overcast 
distributions,  the  results  of  the  cloud  decision  algorithm  will  be  dependent  on  how  the  threshold 
is  applied.  In  this  case,  a  threshold  of  0.18  will  misidentify  some  potentially  cloud-free  stars  as 
being  obscured  by  clouds.  Since  the  number  of  cases  that  fall  into  that  range  appears  small,  the 
results  should  not  be  drastically  affected. 

Algorithm  Results  and  Validation 

When  GetCloudData  is  run,  several  optional  parameters  may  be  included  in  the 
calculation  of  the  LMM  best  fit.  Results  to  date  show  that  there  is  some  dependence  of  the 
results  on  how  these  parameters  are  initialized.  For  instance,  an  array  input  to  MPFit  contains 
limits  for  each  of  the  six  parameters.  Preliminary  studies  showed  that  limiting  the  point-spread 
function  widths  had  an  adverse  affect  on  the  results.  Therefore,  requiring  that  the  background 
and  scaling  factors  be  greater  than  zero  were  the  only  limits  used. 
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Relative  Frequency 


It  was  also  found  that  the  LMM  subroutine  works  best  when  given  initial  parameters  that 
are  close  to  the  minimized  values.  The  background  value  is  set  equal  to  the  minimum  pixel 
value  in  the  input  array.  The  scaling  factor  is  set  equal  to  the  maximum  pixel  value  minus  the 
minimum  pixel  value.  The  x-and  y-widths  are  equal  to  0.5,  since  this  is  a  typical  point-spread 
function  width.  Finally,  the  x-  and  y-centers  are  set  to  the  central  pixel  location. 

An  additional  parameter  that  needs  to  be  investigated  further  is  the  use  of  weighted  pixel 
values.  It  has  been  found  that  leaving  the  pixel  values  unweighted  is  effective  for  low  magnitude 


Figure  2.  Point— Spread  Function  Widths  for  Cloud-Free  Images 
Magnitudes  <  5  (solid)  and  Magnitudes  >  5  (dash) 


Point-Spread  Function  Width  (pixels) 
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Figure  3.  Comparison  of  Star  Contrast  for  Cloud-Free  and  Cloudy  Conditions 


Star  Contrast  under  Cloud-Free  Conditions 
Magnitudes  in  Range  0-5 


Star  Contrast  under  Cloud-Free  Conditions 
Magnitudes  in  Range  5-6 
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Star  Contrast  under  Overcast  Conditions 
Magnitudes  in  Range  0-5 

Star  Contrast  under  Overcast  Conditions 
Magnitudes  in  Range  5-6 

stars,  but  using  weights  that  are  directly  proportional  to  each  pixel  value  is  more  effective  for 
high  magnitude  stars.  The  point-spread  functions  are  more  realistic,  and  the  algorithm  converges 
much  quicker.  One  possible  scenario  is  to  use  unweighted  pixel  values  for  low  magnitude  stars, 
switching  to  weights  beyond  a  certain  magnitude  value.  For  the  data  used  in  this  memo,  no 
weights  were  used. 

Despite  the  large  number  of  dependencies,  LMM  is  quite  robust.  CURVEFIT,  the 
nonlinear  least-squares  routine  included!  with  IDL,  uses  a  gradient  expansion  algorithm.  A 
comparison  of  the  two  methods  using  modeled  stars  in  the  form  of  two-dimensional  Gaussian 
distributions  showed  LMM  converged  to  the  exact  point-spread  function  values  in  most  cases, 
while  the  gradient  expansion  method  did  not.  Results  for  the  two  methods  were  comparable 
when  an  input  array  with  a  width  of  nine  pixels  was  used,  except  for  point-spread  function 
widths  less  than  0.4,  where  LMM  results  were  superior.  Neither  method  worked  well  when  the 
model  star’s  point-spread  function  width  was  much  lower  than  0.3.  When  the  array  width  was 
reduced  to  five  pixels,  the  gradient  expansion  method  performed  far  worse  than  LMM. 

To  validate  the  preliminary  star  algorithm,  several  nighttime  images  were  analyzed. 
Figures  4  to  8  each  show  the  original  image  and  the  image  overlaid  with  a  mask  that  is  color- 
coded  to  represent  the  cases  of  cloud-free  (blue),  overcast  (red),  and  no  usable  stars  (gray).  The 
first  two  examples  are  cloud-free  images  collected  at  0427  UT  on  December  18,  1998  and  0220 
UT  on  December  18,  1997,  respectively.  Figures  6  and  7  are  of  completely  overcast  images 
collected  at  0356  UT  on  November  18,  1998  and  0820  UT  on  January  17,  1999,  respectively. 
Finally,  Figure  8  shows  an  image  partly  covered  with  thin  and  opaque  clouds  collected  at  1127 
UT  on  December  18,  1998. 

Table  1  summarizes  the  number  of  image  cells  out  of  237  that  fall  into  each  category. 
Also  shown  is  the  number  of  image  cells  that  did  not  have  any  stars  with  magnitudes  <5.  While 
the  results  are  not  perfect,  the  overall  strategy  appears  to  be  effective. 


Table  1. 


Number  of  Image  Cells  in  Each  Category 


Image 

Cloud-Free 

Obscured 
by  Cloud 

No  Stars 

#  w/o  Low 
Mag  Stars 

Cloud-Free  1 

228 

8 

1 

38 

Cloud-Free  2 

230 

7 

0 

35 

Overcast  1 

3 

234 

0 

34 

Overcast  2 

1 

235 

1 

41 

Partly-Cloudy 

28 

206 

3 

37 

We  can  also  examine  the  results  by  looking  at  the  stars  individually.  Table  2  gives  the 
algorithm  results  for  catalog  stars  as  a  function  of  magnitude  and  zenith  angle.  The  data  is 
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summarized  from  multiple  cloud-free  images.  The  first  column  lists  the  magnitude  range  for  the 
row.  Results  shown  in  the  other  columns  are  broken  down  by  zenith  angle:  stars  with  zenith 
angles  less  than  70,  50,  and  30°,  respectively.  The  second  column  shows  the  total  number  of 
catalog  stars  in  the  given  magnitude  range.  The  third  column  gives  the  percentage  of  peaks 
identified  as  catalog  stars  using  the  full  algorithm  logic,  including  contrast  and  point-spread 
width  thresholds.  The  last  column  shows  the  mean  residual  in  pixels  between  the  star 
coordinates  found  from  the  Gaussian  best  fit  and  the  image  star  coordinates  corresponding  to  the 
stars  identified  in  column  three. 


Table  2.  Star  Identification  Results  for  Cloud-Free  Scenes  as  a  Function  of 
Star  Magnitude  and  Zenith  Angle 


Number  of  Stars 
(70/50/30°) 

%  Success  by  Zenith  Angle 
(70/50/30°) 

Mean  Pixel  Residual  by 
Zenith  Angle  (70/50/30°) 

0.0-4.0 

296/168  /  75 

0.93  /  0.89  /  0.91 

0.48  /  0.37  /  0.31 

4.0-5.0 

728/421/175 

0.90  /  0.90  /  0.86 

0.51  /  0.43  /  0.35 

5.0-6.0 

2093/1164/438 

0.61  /  0.65  /  0.68 

0.58  /  0.49/0.43 

6.0-7.0 

2104/1151/438 

0.25/0.30/0.30 

0.57  /  0.45/0.43 

It  is  clear  from  the  table  that  the  number  of  available  stars  increases  substantially  with 
magnitude.  It  can  also  be  seen  that  the  ability  to  successfully  identify  stars  decreases  with 
magnitude.  This  is  not  all  that  surprising,  since  the  fainter  stars  are  not  expected  to  stand  out 
from  the  background  as  well  as  brighter  stars,  and  therefore  less  likely  to  be  identified  with  the 
best  fit  algorithm.  The  percentage  of  identified  stars  with  magnitudes  less  than  5,  86-93%,  is 
quite  good.  Beyond  that,  the  percentages  drop  off  substantially,  and  for  magnitudes  6  and 
greater,  the  stars  appear  to  be  a  useless  indicator  for  the  methods  used  here.  Looking  at  the  flags 
that  were  set  for  these  stars,  30%  of  magnitude  5-6  stars  and  56%  of  magnitude  6-7  stars  had 
additional  peaks  nearby.  These  additional  signals  undoubtedly  contributed  to  the  failure  of  the 
best  fit  algorithm  for  most  of  the  cases.  The  residuals  shown  in  the  last  column  increase 
somewhat  with  magnitude,  but  the  change  is  not  significant.  All  cases  had  mean  residuals  less 
than  0.6  pixel. 

While  there  is  a  finite  probability  that,  given  random  pixel  coordinates,  a  peak  associated 
with  a  star  will  be  found  nearby,  examination  of  images  showing  the  location  of  catalog  stars  and 
the  associated  peaks  shows  that  the  method  works  well  for  low  magnitudes.  To  test  the 
likelihood  of  randomly  finding  a  peak  resembling  a  star  in  an  image,  a  routine  called 
GetRandomStars  was  written.  This  routine  uses  a  random  number  generator  to  obtain  dummy 
star  locations  that  are  then  passed  to  the  peak  search  and  best  fit  subroutines,  just  as  with  the 
catalog  stars.  In  the  case  of  completely  cloud-free  images,  there  was  about  20%  chance  that  a 
peak  resembling  a  star  was  located  using  the  star  algorithm,  given  500  random  coordinates. 
Examining  the  resulting  peak  locations  showed  that  the  majority  of  these  peaks  were  indeed 
stars.  The  same  exercise  with  a  completely  overcast  image  resulted  in  less  than  1%  chance  that  a 
peak  resembling  a  star  was  located. 
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Conclusion  and  Future  Work 


This  first  generation  algorithm  appears  to  work  well  for  skies  with  no  moon.  There  are 
several  approaches  to  further  optimization  that  we  feel  would  be  worth  evaluating.  One 
approach  that  has  showed  reasonable  results  involves  the  ratio  between  the  contrast  for  a  star  in  a 
cloud-free  image  and  the  contrast  for  the  same  star  in  an  unknown  image.  This  method  appears 
to  work  well  in  distinguishing  cloud-free,  thin  cloud  and  opaque  cloud  in  the  unknown  image.  In 
this  case,  a  library  of  cloud-free  contrast  values  for  specific  stars  would  need  to  be  created. 
Unfortunately,  light  from  man-made  and  other  sources  creates  some  dependence  of  the  sky 
background  on  zenith  angle.  This  dependence  must  be  removed  before  the  contrast  ratio  method 
could  be  applied.  In  addition,  it  may  be  possible  to  further  optimize  this  algorithm  with 
techniques  such  as  customizing  the  star  list  to  remove  stars  that  are  close  to  a  brighter  star. 
Removing  secondary  peaks  near  the  star  center  from  the  input  array  may  aid  in  calculating  the 
best  fit  Gaussian,  particularly  with  high  magnitude  stars.  The  algorithm  also  needs  to  be 
evaluated  using  data  from  other  sites  and  using  raw  data  as  opposed  to  calibrated  data. 

However,  we  feel  that  this  first  generation  algorithm  is  sufficiently  accurate  to  merit 
converting  to  C-code  that  can  be  included  in  field  software.  Once  a  usable  field  version  is 
written,  future  developments  could  include  determination  of  optical  depth  from  the  calibrated 
images,  evaluating  the  transition  to  partial  moonlight,  and  increasing  the  angular  extent  and/or 
angular  resolution  of  the  algorithm. 
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Grey 


Figure  5.  Cloud-Free  Scene 
Red  =  >50%  Failure,  Blue  =  <50%  Failure 
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Figure  7.  Overcast  Scene 
Grey  =  No  Stars,  Red  =  >50%  Failure,  Blue 


<50%  Failure 
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Figure  8.  Partly-Cloudy  Scene 

Grey  =  No  Stars,  Red  =  >50%  Failure,  Blue  =  <50%  Failure 
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